home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / earthmap / _earthmap / maps / program / c / domap
Encoding:
Text File  |  1989-04-01  |  11.9 KB  |  566 lines
  1. /*
  2.  * Usage:
  3.  * domap [file1] [file2] [file3] [...] [options] | pen xcenter=0 ycenter=0
  4.  *                           ^^^
  5.  *                           A Vplot filter
  6.  *
  7.  * If project=round, this program projects the Earth onto a plane tangent
  8.  * to its surface, as seen from your "eye", where the position of your eye
  9.  * is given by the following parameters:
  10.  * 
  11.  * If project=round
  12.  *    lat=37.45       <- Latitude of eye
  13.  *    long=-122.21    <- Longitude of eye
  14.  *    dist=100.       <- Distance of eye above Earth, in miles
  15.  *    inch=110.          <- Miles per inch at center of projection
  16.  *
  17.  * If project=flat, a crude latitude-longitude plot is made instead.
  18.  *
  19.  * If project=flat
  20.  *    long=-122.21    <- Longitude centered in the plot
  21.  *
  22.  * file1, file2, file3, etc, are map information files
  23.  * such as are in ./Namer.Map/*.cbd, etc.
  24.  *
  25.  * A slew of other options are used to control what sorts of objects
  26.  * get plotted, and in what colors. These are:
  27.  * 
  28.  * pby1=6   <- Vplot color to make this type of object; 0 means omit
  29.  * pby2=6
  30.  * pby3=5
  31.  * 
  32.  * bdy1=7
  33.  * bdy2=7
  34.  * bdy3=6
  35.  * 
  36.  * riv1=5
  37.  * riv2=1
  38.  * riv3=1
  39.  * riv4=1
  40.  * riv5=5
  41.  * riv6=1
  42.  * riv7=1
  43.  * riv8=1
  44.  * riv9=1
  45.  * riv10=1
  46.  * riv11=1
  47.  * riv12=1
  48.  * 
  49.  * cil1=7
  50.  * cil2=4
  51.  * cil3=1
  52.  * cil4=1
  53.  * cil5=1
  54.  * cil6=1
  55.  * cil7=1
  56.  * cil8=1
  57.  * cil9=1
  58.  * cil10=1
  59.  * cil11=1
  60.  * cil12=1
  61.  * cil13=1
  62.  * cil14=1
  63.  * cil15=1
  64.  *
  65.  * The meanings of these different objects are defined in the
  66.  * "Database_info" file.
  67.  *
  68.  * To read the data files properly, you may need to rewrite the
  69.  * "twiddle32" and "twiddle16" subroutines. The data files are
  70.  * integers, and so preoper reading of them depends on the byte
  71.  * ordering of your machine.
  72.  *
  73.  */
  74. #include <stdio.h>
  75. #include <math.h>
  76. #include <vplot.h>
  77. #include <strings.h>
  78.  
  79. /*
  80.  * Much of this program was written by Brian Reid at the Dec
  81.  * Western Research Lab. The rest was written by Joe Dellinger
  82.  * at Stanford University. Neither of us much cares what becomes
  83.  * of it. Sat Apr  1 01:26:51 PST 1989
  84.  */
  85.  
  86. /* ----------------------------------------------------------------------
  87.  * Read and process one binary map file.
  88.  */
  89.  
  90. #include <sys/file.h>
  91.  
  92. #define SHORTFLAG 0x4000
  93. #define CBD_MAGIC 0x20770002
  94. #define BIT32 int
  95. #define BIT16 short
  96.  
  97. #define CONV (3.141592654 / 180.)
  98.  
  99. /*
  100.  * Input parameters:
  101.  * viewlat, viewlong = lat, long in center of picture;
  102.  * viewdist = distance of eye from earth, in miles;
  103.  * viewscale = miles per inch at center of picture
  104.  *
  105.  * For project=round, all of these are ignored except
  106.  * viewlong.
  107.  */
  108. double          viewlat, viewlong, viewdist, viewscale;
  109. char            projection[80];
  110.  
  111. /*
  112.  * for getpar
  113.  */
  114. int             xargc;
  115. char          **xargv;
  116.  
  117. struct cbdhead
  118. {
  119.     BIT32           magic;    /* Magic number */
  120.     BIT32           dictaddr;    /* Offset of segment dictionary in file */
  121.     BIT32           segcount;    /* Number of segments in file */
  122.     BIT32           segsize;    /* Size of segment dictionary */
  123.     BIT32           segmax;    /* Size of largest segment's strokes, /2 */
  124.     BIT32           fill[(40 / sizeof (BIT32)) - 5];    /* Filler */
  125. };
  126.  
  127. main (argc, argv)
  128.     int             argc;
  129.     char          **argv;
  130. {
  131. int             Ifile, segcount, idx, idy, segbufsize, olt, oln, j, k, jstroke, iseg;
  132. int             colfile[20];
  133. char            string[80], string2[80], *sptr;
  134. char            filename[20][80];
  135. int             filenum;
  136. BIT32           i32;
  137. BIT16          *segbuf;
  138. double          x, y, z, x1, x2, xc, y1, y2, yc, z1, z2, zc;
  139. double          dist, cent, clipdist, horizon, my_acos ();
  140. double          longi, lati;
  141.  
  142. struct segdict
  143. {
  144.     BIT32           segid;
  145.     BIT32           maxlat, minlat, maxlong, minlong;
  146.     BIT32           absaddr;
  147.     BIT16           nbytes;
  148.     BIT16           rank;
  149. }              *sd, *sdbuf;
  150.  
  151. struct segment
  152. {
  153.     BIT32           orgx, orgy;
  154.     BIT32           id;
  155.     BIT16           nstrokes;
  156.     BIT16           dummy;
  157. }               sb;
  158.  
  159. struct cbdhead  header;
  160.  
  161. #ifdef DEBUG
  162. int             deg, min, sec, count;
  163. #endif
  164.  
  165.  
  166.     xargc = argc;
  167.     xargv = argv;
  168.  
  169.     vp_filep (stdout);
  170.  
  171.     viewlat = 40.;
  172.     viewlong = -90.;
  173.     viewdist = 100000.;
  174.     viewscale = 100.;
  175.  
  176.     getpar ("lat", "g", &viewlat);
  177.     getpar ("long", "g", &viewlong);
  178.     getpar ("dist", "g", &viewdist);
  179.     getpar ("inch", "g", &viewscale);
  180.  
  181.     viewlat *= CONV;
  182.     viewlong *= CONV;
  183.     viewdist /= 3950.;
  184.     viewscale /= 3950.;
  185.  
  186.     strcpy (projection, "round");
  187.     getpar ("project", "s", projection);
  188.  
  189.     if (projection[0] == 'r')
  190.     vp_scale (1. / viewscale, 1. / viewscale);
  191.     else
  192.     vp_scale (5. / (180. * CONV), 5. / (180. * CONV));
  193.  
  194.     clipdist = my_acos (1. / (1. + viewdist));
  195.  
  196.     vp_color (WHITE);
  197.     vp_penup ();
  198.     vp_fat (1);
  199.  
  200.     if (projection[0] == 'r')
  201.     {
  202.     horizon = tan (clipdist / 2.);
  203.     for (j = 0; j <= 360 * 10; j++)
  204.     {
  205.         vp_upendn (horizon * cos (j * CONV / 10.), horizon * sin (j * CONV / 10.));
  206.     }
  207.     }
  208.     else
  209.     {
  210.     vp_move (-5., -2.5);
  211.     vp_draw (5., -2.5);
  212.     vp_draw (5., 2.5);
  213.     vp_draw (-5., 2.5);
  214.     vp_draw (-5., -2.5);
  215.     }
  216.  
  217.     for (j = 0; j < 20; j++)
  218.     {
  219.     strcpy (filename[j], "");
  220.     }
  221.  
  222.     for (filenum = 1; filenum < argc; filenum++)
  223.     {
  224.  
  225.     strcpy (filename[filenum], argv[filenum]);
  226.  
  227.     if (filename[filenum][0] == '\0'
  228.         ||
  229.         index (filename[filenum], '=') != NULL
  230.      )
  231.         continue;
  232.  
  233.     Ifile = open (filename[filenum], O_RDONLY, 0);
  234.  
  235.     if (Ifile == -1)
  236.         continue;
  237.  
  238.     for (j = 0; j < 20; j++)
  239.     {
  240.         colfile[j] = 0;
  241.     }
  242.  
  243.     sptr = rindex (filename[filenum], '.');
  244.     if (sptr == NULL)
  245.         continue;
  246.     if (strcmp (sptr, ".cbd") != 0)
  247.         continue;
  248.  
  249.     strcpy (string2, sptr - 3);
  250.     string2[3] = '\0';
  251.  
  252.     colfile[1] = WHITE;
  253.     colfile[2] = CYAN;
  254.  
  255.     for (j = 0; j < 20; j++)
  256.     {
  257.         sprintf (string, "%s%d", string2, j);
  258.         getpar (string, "d", &colfile[j]);
  259.     }
  260.  
  261. /*
  262.  * Read in the file header, check it for the correct magic number,
  263.  * and learn the address of the segment dictionary
  264.  */
  265.     read (Ifile, &header, sizeof header);
  266.  
  267.     twiddle32 (&(header.magic));
  268.     twiddle32 (&(header.dictaddr));
  269.     twiddle32 (&(header.segcount));
  270.     twiddle32 (&(header.segsize));
  271.     twiddle32 (&(header.segmax));
  272.  
  273.     if (header.magic != CBD_MAGIC)
  274.     {
  275.         fprintf (stderr, "File has bad magic number %d != %d\n",
  276.              header.magic, CBD_MAGIC);
  277.         exit (2);
  278.     }
  279.  
  280. /* allocate space for the segment buffer */
  281.     segbufsize = 2 * header.segmax;
  282.     segbuf = (BIT16 *) malloc (50 + segbufsize);
  283.  
  284. /* allocate space for the segment dictionary */
  285.     sdbuf = (struct segdict *) malloc (100 + header.segsize);
  286.     sd = sdbuf;
  287.     sd++;
  288.  
  289. /* Go read in the segment dictionary (it's at the end of the file) */
  290.     lseek (Ifile, header.dictaddr, L_SET);
  291.     j = read (Ifile, sd, header.segsize);
  292.  
  293.     if (j < header.segsize)
  294.     {
  295.         fprintf (stderr, "File segment dictionary expecting %d got %d\n",
  296.              header.segsize, j);
  297.         close (Ifile);
  298.         free (sdbuf);
  299.         free (segbuf);
  300.         return;
  301.     }
  302.  
  303. /*
  304.  * Now look at each segment and decide if we're
  305.  * going to keep it or not
  306.  */
  307.     segcount = header.segcount;
  308.  
  309.     sd = sdbuf;
  310.     for (iseg = 1; iseg <= segcount; iseg++)
  311.     {
  312.         sd++;        /* does this really work? wow! */
  313.  
  314.         twiddle32 (&(sd->segid));
  315.         twiddle32 (&(sd->maxlat));
  316.         twiddle32 (&(sd->minlat));
  317.         twiddle32 (&(sd->maxlong));
  318.         twiddle32 (&(sd->minlong));
  319.         twiddle32 (&(sd->absaddr));
  320.         twiddle16 (&(sd->nbytes));
  321.         twiddle16 (&(sd->rank));
  322.  
  323.         if (colfile[sd->rank] == 0)
  324.         continue;
  325.  
  326.         if (projection[0] == 'r')
  327.         {
  328.         longi = sd->maxlong * CONV / 3600.;
  329.         lati = sd->maxlat * CONV / 3600;
  330.  
  331.         x = sin (longi - viewlong) * cos (lati);
  332.         y = cos (longi - viewlong) * cos (lati);
  333.         z = sin (lati);
  334.  
  335.         y1 = y * cos (viewlat) + z * sin (viewlat);
  336.         z1 = -y * sin (viewlat) + z * cos (viewlat);
  337.         x1 = x;
  338.  
  339.         longi = sd->minlong * CONV / 3600;
  340.         lati = sd->minlat * CONV / 3600;
  341.         x = sin (longi - viewlong) * cos (lati);
  342.         y = cos (longi - viewlong) * cos (lati);
  343.         z = sin (lati);
  344.  
  345.         y2 = y * cos (viewlat) + z * sin (viewlat);
  346.         z2 = -y * sin (viewlat) + z * cos (viewlat);
  347.         x2 = x;
  348.  
  349.         xc = (x1 + x2) / 2.;
  350.         yc = (y1 + y2) / 2.;
  351.         zc = (z1 + z2) / 2.;
  352.         x = sqrt (xc * xc + yc * yc + zc * zc);
  353.         xc /= x;
  354.         yc /= x;
  355.         zc /= x;
  356.  
  357.         dist = my_acos (xc * x1 + yc * y1 + zc * z1);
  358.         cent = my_acos (yc);
  359.  
  360.         if (cent - dist > clipdist)
  361.             continue;
  362.         }
  363.  
  364.         vp_color (colfile[sd->rank]);
  365.         if (sd->rank == 1)
  366.         vp_fat (1);
  367.         else
  368.         vp_fat (0);
  369.  
  370.         lseek (Ifile, sd->absaddr, L_SET);
  371.         read (Ifile, &sb, sizeof sb);
  372.         twiddle32 (&(sb.orgx));
  373.         twiddle32 (&(sb.orgy));
  374.         twiddle32 (&(sb.id));
  375.         twiddle16 (&(sb.nstrokes));
  376.  
  377.         if (sd->nbytes > segbufsize)
  378.         {
  379.         fprintf (stderr, "Segment %d needs %d bytes; buffer limit is %d.\n", iseg, sd->nbytes, segbufsize);
  380.         close (Ifile);
  381.         free (sdbuf);
  382.         free (segbuf);
  383.         return;
  384.         }
  385.         read (Ifile, segbuf, sd->nbytes);
  386.  
  387.         k = 0;
  388.  
  389.         oln = sb.orgx;
  390.         olt = sb.orgy;
  391.  
  392. #ifdef DEBUG
  393.         count = 1;
  394.         deg = abs (oln) / 3600;
  395.         min = (abs (oln) - deg * 3600) / 60;
  396.         sec = abs (oln) - deg * 3600 - min * 60;
  397.         if (oln > 0)
  398.         fprintf (stderr, "+ %d %d %d   %d\n", deg, min, sec, count);
  399.         else
  400.         fprintf (stderr, "- %d %d %d   %d\n", deg, min, sec, count);
  401. #endif
  402.  
  403.         vp_penup ();
  404.         project (CONV * olt / 3600., CONV * oln / 3600.);
  405.  
  406.         for (jstroke = 1; jstroke <= sb.nstrokes; jstroke++)
  407.         {
  408.         twiddle16 (&segbuf[k]);
  409.         if (segbuf[k] & SHORTFLAG)
  410.         {
  411. #ifdef DEBUG
  412.             fprintf (stderr, "s ");
  413. #endif
  414. /* Flag bit on: unpack a 16-bit field into dx and dy */
  415.             i32 = segbuf[k++];
  416.             if (i32 > 0)
  417.             i32 &= ~SHORTFLAG;
  418.             idy = i32 & 0xFF;
  419.             if (idy & 0x80)
  420.             idy |= ~0xFF;    /* extend sign */
  421.             idx = i32 >> 8;
  422.             if (idx & 0x80)
  423.             idx |= ~0xBF;    /* extend sign */
  424.         }
  425.         else
  426.         {
  427. #ifdef DEBUG
  428.             fprintf (stderr, "L ");
  429. #endif
  430. /* Flag bit off: take dx and dy from 32-bit fields. */
  431.             idx = segbuf[k++];
  432.             if (idx < 0)
  433.             idx |= SHORTFLAG;
  434.             twiddle16 (&segbuf[k]);
  435.             idx = (idx << 16) | segbuf[k];
  436.             k++;
  437.  
  438.             twiddle16 (&segbuf[k]);
  439.             idy = segbuf[k];
  440.             k++;
  441.             if (idy < 0)
  442.             idy |= SHORTFLAG;
  443.             twiddle16 (&segbuf[k]);
  444.             idy = (idy << 16) | segbuf[k];
  445.             k++;
  446.         }
  447.         olt = olt + idy;
  448.         oln = oln + idx;
  449.  
  450. #ifdef DEBUG
  451.         count++;
  452.         deg = abs (oln) / 3600;
  453.         min = (abs (oln) - deg * 3600) / 60;
  454.         sec = abs (oln) - deg * 3600 - min * 60;
  455.         if (oln > 0)
  456.             fprintf (stderr, "  + %d %d %d   %d\n", deg, min, sec, count);
  457.         else
  458.             fprintf (stderr, "  - %d %d %d   %d\n", deg, min, sec, count);
  459. #endif
  460.  
  461.         project (CONV * olt / 3600., CONV * oln / 3600.);
  462.         }
  463.     }
  464.     close (Ifile);
  465.     free (sdbuf);
  466.     free (segbuf);
  467.     }
  468. }
  469.  
  470. /*
  471.  * Change the two routines below to accomodate the byte ordering on your
  472.  * machine! - Joe Dellinger
  473.  */
  474.  
  475. twiddle32 (i32)
  476.     BIT32          *i32;
  477. {
  478. BIT32           temp;
  479.  
  480.     temp = 0;
  481.     temp |= 0xFF000000 & ((0x000000FF & *i32) << 24);
  482.     temp |= 0x00FF0000 & ((0x0000FF00 & *i32) << 8);
  483.     temp |= 0x0000FF00 & ((0x00FF0000 & *i32) >> 8);
  484.     temp |= 0x000000FF & ((0xFF000000 & *i32) >> 24);
  485.     *i32 = temp;
  486. }
  487.  
  488. twiddle16 (i16)
  489.     BIT16          *i16;
  490. {
  491. BIT16           temp;
  492.  
  493.     temp = 0;
  494.     temp |= 0xFF00 & ((0x00FF & *i16) << 8);
  495.     temp |= 0x00FF & ((0xFF00 & *i16) >> 8);
  496.     *i16 = temp;
  497. }
  498.  
  499.  
  500. project (lati, longi)
  501.     double          lati, longi;
  502. {
  503. double          xx, yy, zz;
  504. double          x, y, z;
  505. double          dot, factor;
  506.  
  507. static double   xlast = 0.;
  508.  
  509.     if (projection[0] == 'r')
  510.     {
  511.     x = sin (longi - viewlong) * cos (lati);
  512.     y = cos (longi - viewlong) * cos (lati);
  513.     z = sin (lati);
  514.  
  515.     yy = y * cos (viewlat) + z * sin (viewlat);
  516.     zz = -y * sin (viewlat) + z * cos (viewlat);
  517.     xx = x;
  518.  
  519.     dot = xx * xx + zz * zz - yy * (viewdist + (1. - yy));
  520.  
  521.     if (dot > 0.)
  522.     {
  523.         vp_penup ();
  524.         return;
  525.     }
  526.  
  527.     factor = viewdist / (viewdist + (1. - yy));
  528.     xx *= factor;
  529.     zz *= factor;
  530.  
  531.     vp_upendn (xx, zz);
  532.     }
  533.     else
  534.     {
  535.     x = longi - viewlong;
  536.     while (x >= 180 * CONV)
  537.         x -= 360 * CONV;
  538.     while (x < -180 * CONV)
  539.         x += 360 * CONV;
  540.     y = lati;
  541.  
  542.     if (fabs (xlast - x) > 180. * CONV)
  543.         vp_penup ();
  544.     xlast = x;
  545.  
  546.     vp_upendn (x, y);
  547.     }
  548. }
  549.  
  550. double
  551. my_acos (x)
  552.     double          x;
  553. {
  554.     if (fabs (x) >= 1.)
  555.     return 0.;
  556.     else
  557.     return acos (x);
  558. }
  559.  
  560. /*
  561.  * Just eat any error message from getpar.
  562.  */
  563. err ()
  564. {
  565. }
  566.